clc; clearvars -except NumberConsecutiveshocks LogLevel shockchoice
   LogLevel  =  1;
%  frequency = 12; % Monthly
% %frequency =  4; % Quarterly
%-------------------------------------------------------------------------%
    load RVSP500.mat
   %sample=[ 1983 2; 2014  12];%as in  Berger, Dew-Becker and Giglio (2017), "Uncertainty shocks as second moment
   %
    sample=[ 1983 2; 2019  12];
   %sample=[ 1990 1; 2019   6]; %as in Sahalia et al. (2020) ``Volatility and Uncertainty Disconnected''
   %sample=[ 2007 7; 2009   6]; %as in Sahalia et al. (2020) ``Volatility and Uncertainty Disconnected''

% take subset
startsample=find(yyyyRV==sample(1,1)&mmRV==sample(1,2));
endsample  =find(yyyyRV==sample(2,1)&mmRV==sample(2,2));

dates_RV = yyyyRV(startsample:endsample)*100+mmRV(startsample:endsample);
RV = rv(startsample:endsample);
% Replicating Sahalia et al. (2020) ``Volatility and Uncertainty Disconnected''
fprintf('Realized volatility Market Returns %3.3f\n',mean(100*sqrt(RV*12)))

%% CHECK -- From BDBG ReStud
% load('alldata.mat','s_cme_rvar','s_cme_vixF','cme_RV','masterdates','masterdates_matlab')
% f=find(floor(masterdates/100)>=1983);
%
% load('RVSP500.mat')
% rvSP=rv;
% clear rv
% load('RVMKT.mat')
% cstr=cal(1926,7,12,size(rv));
% ss = ical(1983,1,cstr); se=ical(2014,12,cstr);
% plot([         cme_RV(f)       rvSP(ss:se)])
% plot([    sqrt(cme_RV(f)) sqrt(rvSP(ss:se))]*sqrt(12)*100)
% plot([                    sqrt(rvSP(ss:se))*sqrt(12)*100 s_cme_rvar(f)])
% plot([    sqrt(cme_RV(f))*sqrt(12)*100                   s_cme_rvar(f)])
%
% plot([log(sqrt(cme_RV)*sqrt(12)*100) ls_cme_rvar])

if LogLevel
    ls_RV = log(sqrt(RV)*sqrt(12)*100);
else
     s_RV =     sqrt(RV)*sqrt(12)*100 ;
end

clearvars -except sample NumberConsecutiveshocks RV ls_RV SP500 LogLevel dates_RV sample shockchoice

dates_RV = dates_RV(2:end);
%% Load VIX
 VIX = xlsread('VIXCLS.xls','VIXmerged');
 VIXdates = VIX(1:end,1);
 VIX      = VIX(1:end,2); % plot([exp(ls_RV) VIX])
startsample=find(VIXdates==(sample(1,1)*100+sample(1,2)));
endsample  =find(VIXdates==(sample(2,1)*100+sample(2,2)));
 if LogLevel
     logVIX = log(VIX(startsample:endsample));
 else
     logVIX =     VIX(startsample:endsample);
 end
 clear startsample endsample VIXdates

 %% Set-up VAR as in BdBG ReStud
sample=[ 1983 2; 2019  12]; 
 data = xlsread('../../Dataset/VARdataBBupdatedMonthly.xlsx','MainData');
startsample=find(data(:,1)==sample(1,1)*100+sample(1,2));
endsample  =find(data(:,1)==sample(2,1)*100+sample(2,2));

data = data(startsample:endsample,2:end);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATA TRANSFORMATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data(:,[3:11]) = log(data(:,[3:11]))*100;
%--------------------------------------%

%All VARs below use four lags, as suggested by the Akaike information criterion for the main specification.
%The vector of variables in the benchmark specificationis rv,v,FFR,ip,emp,

 dataMatrix= [ls_RV  logVIX   data(:,1)                  data(:,3)                data(:,10)]; %Total Nonfarm Payroll (PAYEMS)
%dataMatrix= [ls_RV  logVIX   data(:,1)                  data(:,3)                data(:,11)]; %All Employees, Manufacturing (MANEMP)
dataName   = {  'rv',  'VIX','Short-term nominal rate', 'Industrial production', 'Employment'}; 
clear data
VarsToPlot = [4 5]; pos_shock = 1;

nobs       = size(dataMatrix,1);

%-------------------------------------------------------------
% Step 0: parameter settings
%-------------------------------------------------------------
%nlags  = 3;                   % number of lags in VAR 
 nlags  = 4;                   % number of lags in VAR 
 nconst = 1;                   % 3 = for constant, linear and quadratic trend, 2 = for constant and trend, 1 = for constant, 0 otherwise 

%-------------------------------------------------------------------------%
prior=1;                      % for prior,  1   = standard prior
                              %             inf =     flat prior

%set prior for first lag of each variable (if estimation with Minnesota prior)
[~,coly]=size(dataMatrix);
lev=ones(coly,1);                                          
%-------------------------------------------------------------------------%
% nimp = 12;                  % horizon of VD and IRFs to be shown         
  nimp = 36;                  % horizon of VD and IRFs to be shown         

est_method = 0;               % 0 : OLS point estimation
                              % 1 : Bayesian estimation with Minnesota prior 
                              % 2 : OLS point estimation and classical bootstrap    
                              
%-----------------------------------------------------------------------
% Step 1: load and transform data
%-----------------------------------------------------------------------
[y,x] = mk_vdm(dataMatrix,nlags,0);
%y is T x nvars matrix of rhs variables
%x is T x nvars*nlags  of lhs variables
[T,nvars]=size(y);

%---------------------------------------------------------------------
% Step 3: Estimate VAR, extract shocks, FEV shares explained, and IRFs 
%---------------------------------------------------------------------
if nconst==1
   const = ones(T,1); 
   x = [x const];
   ncoeffs = nvars*nlags+1;
elseif nconst==2
   x = [x ones(T,1) (1:T)'];
   ncoeffs = nvars*nlags+2;
elseif nconst==3
    x = [x ones(T,1) (1:T)' ((1:T).^2)'];
%   dummy_1975Q2 = dates_q_==1975.5;
%   x = [x ones(T,1) (1:T)' ((1:T).^2)' dummy_1975Q2(1:end-nlags)];
   ncoeffs = nvars*nlags+3;   
else
   ncoeffs = nvars*nlags; 
end

%ncoeffs + 1 x nvars matrix of regression coefficients
b=x\y;       %OLS coefficients: ncoeffs x nvars (** more efficient way of computing least squares; i.e. b=inv(x'*x)*x'*y;
res = y - x*b;          %T x nvar matrix of residuals
vmat = (1/T)*(res'*res);
stds=sqrt(diag(vmat));

%OLS point estimates
[FEVD_IRF,v]  =RecursivenessIdentification(b,res,vmat,nvars,nlags,nimp, pos_shock);

%----------------------------------------------------------------------
% Step 3: Plots
%----------------------------------------------------------------------
% clearvars -except dataName VarsToPlot FEVD_IRF plottype nvars
% save(strcat('IRFDATA_',plottype))

% %FEV shares explained by shock
% plotFEV(FEVD_IRF,vars,slope,use_yields,shock_extract);
% 
%IRFs to shock
% plotIRF(FEVD_IRF,vars,slope,use_yields,shock_extract)
REPLICATION_BdBG = 0;
if REPLICATION_BdBG
% FEVD_IRF
[nimp,~,~]=size(FEVD_IRF);
% R=3; %round(N/2);
h=figure('Color',[0.9412 0.9412 0.9412],'Position',[1 1 800-100 600-100],'Name','IRFs to 1st shock');
figure(h);
plotnum=0;
for n=VarsToPlot
    plotnum=1+plotnum;
    if     length(VarsToPlot)==2 
        subplot(1,2,plotnum);             
    elseif length(VarsToPlot)==5
        subplot(2,3,plotnum);   
    end
    if ndims(FEVD_IRF)==3   %if available, plot confidence bands
        grpyat = [(1:nimp)' FEVD_IRF(1:nimp,nvars+n,2); (nimp:-1:1)' FEVD_IRF(nimp:-1:1,nvars+n,3)];
        patch(grpyat(:,1),grpyat(:,2),[0.7 0.7 0.7],'edgecolor',[0.65 0.65 0.65]);
    end
%     end
    hold on;
% /// responses of industrial production and employment to a 90-point upward EPU innovation,
% ~ 4 std dev
%   plot(1:nimp,       FEVD_IRF(:,nvars+n,1)                      ,'k-');%,'LineWidth',2);
    plot(1:nimp,3.7292*FEVD_IRF(:,nvars+n,1)                      ,'k-');%,'LineWidth',2);
%     end
%   plot(1:nimp,FEVD_IRF(:,nvars+n,1)/FEVD_IRF(1,nvars+2,1),'k-');%,'LineWidth',2);
    plot(1:nimp,zeros(nimp),':k');
    title(dataName(n),'FontSize',14) %'FontWeight','bold',
    ylabel('percent' ,'FontSize',12)
    xlabel('Months','FontSize',12)
    axis tight
   %axis([0 nimp ylow(n) yhigh(n) ]);
    hold off;
end
end
uRV = v(1:end);
ls_RV    =    ls_RV(nlags+1:end);
logVIX   =   logVIX(nlags+1:end);
clearvars -except NumberConsecutiveshocks ls_RV uRV FromM2Q logVIX LogLevel shockchoice

%% FIT AR(1)  --> Shocks 1983/M6 -- 2019/M12    
%    NumberConsecutiveshocks = 2;
% %  NumberConsecutiveshocks = 4;
% /////////////////////////////////////////////////////////////////////////
%----------- Construct indicators for seq. of shocks with 2 ident schemes-%
IndPreviousNShocksRV= zeros(size(uRV));
% /////////////////////////////////////////////////////////////////////////
for tt=NumberConsecutiveshocks:length(uRV)-1 
    if tt==NumberConsecutiveshocks
%       if (                                      sum(uEPU(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPU(tt+1)<0);  IndPreviousNShocksEPU(tt) = 1; end
        if (                                      sum(uRV(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                );  IndPreviousNShocksRV(tt) = 1; end
    else
%       if (uRV(tt-NumberConsecutiveshocks)<0 && sum(uRV(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uRV(tt+1)<0);   IndPreviousNShocksRV(tt) = 1; end
        if (uRV(tt-NumberConsecutiveshocks)<0 && sum(uRV(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks               );   IndPreviousNShocksRV(tt) = 1; end
    end
end
clear tt
% %% To quarterly
% horizonAggreg = 3;
%   logEPUQ = logRV(13:horizonAggreg:end);
%   
%   uEPUQ = chgfreq(uRV,horizonAggreg,horizonAggreg,0); uEPUQ = uEPUQ/std(uEPUQ);
% 
% % /////////////////////////////////////////////////////////////////////////
% %----------- Construct indicators for seq. of shocks with 2 ident schemes-%
% IndPreviousNShocksRVQ= zeros(size(uEPUQ));
% % /////////////////////////////////////////////////////////////////////////
% for tt=NumberConsecutiveshocks:length(uEPUQ)-1 
%     if tt==NumberConsecutiveshocks
% %       if (                                       sum(uEPUQ(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPUQ(tt+1)<0);  IndPreviousNShocksRV(tt) = 1; end
%         if (                                       sum(uEPUQ(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                 );  IndPreviousNShocksRVQ(tt) = 1; end
%     else
% %       if (uEPUQ(tt-NumberConsecutiveshocks)<0 && sum(uEPUQ(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks && uEPUQ(tt+1)<0);  IndPreviousNShocksRV(tt) = 1; end
%         if (uEPUQ(tt-NumberConsecutiveshocks)<0 && sum(uEPUQ(tt-NumberConsecutiveshocks+1:tt)>0)==NumberConsecutiveshocks                 );  IndPreviousNShocksRVQ(tt) = 1; end
%     end
% end
% clear tt
%